Dynamic transition in an atomic glass former: a molecular dynamics evidence 
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We find that a Lennard- Jones mixture displays a dynamic phase transition between an active 
regime and an inactive one. By means of molecular dynamics simulations and of a finite-size study, 
we show that the space time dynamics in the supercooled regime coincides with a dynamic first 
order transition point. 



In order to realize that a physical system has fallen 
into a glassy state one must either drive the system out 
of equilibrium, or investigate its relaxation properties. 
It is indeed a well-known fact that since no static sig- 
nature is available, one has to resort to nonequilibrium 
protocols or measurements to identify a glassy state. The 
supercooled regime, which sets in before the material ac- 
tually freezes into a solid glass, is also characterized by 
anomalous temporal properties, such as the increase in 
the viscosity, and ageing phenomena during relaxation 
processes. The density-density autocorrelation function 
at a microscopic scale, instead of exhibiting a straightfor- 
ward exponential relaxation, develops a plateau (over a 
duration conventionally called r a ) before relaxing to zero. 
The present work is concerned with the dynamical prop- 
erties in the supercooled regime. The idea that dynam- 
ical heterogeneities, long-lived large scale spatial struc- 
tures, are the trademark of the slow and intermittent dy- 
namics in the supercooled regime dates back to almost 
thirty years [l[ and many developments have occurred 
since Here, instead of resorting to local probes such 

as the van Hove function, the self-intermediate scattering 
function, or the nonlinear dynamical susceptibility, - i.e 
on multi-point correlation functions in space and time-, 
we prefer to rely on a global characterization of the sys- 
tem's dynamics. 

More recently, it was advocated by Garrahan et aZ.Q 
that, at the level of Kinctically Constrained Models 
(KCM), dynamical heterogeneities appeared as the con- 
sequence of a universal phase transition mechanism, 
largely independent of the specific model under consid- 
eration. The phase transitions at work are not of a 
conventional type: they occur in the space of trajecto- 
ries of realizations the system has followed over a large 
time duration, instead of occurring in the space of avail- 
able configurations, as the liquid-gas transition, or a 
putative thermodynamic glass transition, for instance, 
do. There exists a well-established body of mathemati- 
cal physics literature to analyze these phase transitions, 
which is based, in spirit, on the thermodynamic formal- 
ism of Ruelle ja, [?J , and which was adapted [|[ and then 
exploited [f| in numerical and analytical studies of the 



KCM's. The idea behind these studies is to follow the 
statistics, upon sampling the physical realizations of the 
system over a large time duration, of a space and time 
extensive observable. This observable characterizes the 
level of dynamical activity in the course of the system's 
time evolution. At this stage, the term activity refers 
to the intuitive picture that an active trajectory is char- 
acterized by many collective rearrangements needed to 
escape from energy basins, whereas an inactive trajec- 
tory is dominated by the rapid in-cage rattling motion 
of the particles without collective rearrangements. This 
activity observable is used to partition trajectories into 
two groups -the inactive and the active ones- depending 
on whether their activity is above or below the average 
activity. 

An important step forward in probing the relevance 
of the dynamic phase transition scenario to realistic 
glasses is the work of Hedges et al. in which an 
atomistic model was considered. The system studied 
by these authors is a mixture of two species of particles 
interacting via a Lennard- Jones potential otherwise 
known as the Kob- Andersen (KA) model [13], and 
which has been shown to fall easily into a glassy state 
without crystallizing. Endowed with a Monte-Carlo or 
a Newtonian dynamics, Hedges et al. implemented the 
Transition Path Sampling (TPS) method to produce the 
probability distribution function (pdf) of the activity, at 
finite times. These authors succeeded in showing that for 
a finite system, as the observation time was increased, 
the pdf of the activity develops a bimodal structure. In 
the light of previous works on lattice models for glass 
formers, they interpreted their result as an evidence 
for a bona fide dynamic phase transition, expected to 
occur once the infinite size and infinite observation time 
limits are reached. While the work of Hedges et al. is a 
definite breakthrough towards a better understanding of 
atomistic glass formers, it leaves a number of questions 
unanswered. These are the subject of the present work. 
Dynamic phase transitions occur in the large time limit, 
as exemplified by the pioneering works of Ott et al. fll| , 
Cvitanovic fl2l or Szepfalusy and Tel |l3[ (see Beck 
and Schlogl [14J for further references). The dynamic 
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transitions at work in glass formers are of a subtler form, 
since they only emerge upon, in addition, considering 
the infinite system size limit. 

Our first goal in this letter is to show the existence of 
a phase transition. Our second goal is to be more precise 
about the location of the phase transition itself. This is 
indeed a central issue: if the critical point is away from 
the typical measurable activity then the transition sce- 
nario is at best a crossover. If, on the contrary, it is shown 
that the phase transition actually takes place for typical 
trajectories then experimentally observable consequences 
are expected. Our third goal is to provide a rigorous 
definition of the activity. It is our feeling that deeper 
insight will be gained if the notion of activity can be 
characterized by standard physical concepts (like forces 
between particles) rather than on phcnomenological con- 
siderations. 

In the present work, we consider a KA mixture, which 
we study by extending the methods of Molecular Dy- 
namics (MD) to the study of temporal large deviations. 
By implementing a MD version of the cloning algorithm 
of Giardina, Kurchan and Peliti [lj| we are able to fol- 
low the statistics of the large deviations of the activity. 
Our activity observable is a physically transparent quan- 
tity which can be related to the rate at which the sys- 
tem escapes from a given location in phase space. As 
we shall soon see, the typical trajectories of the system 
in phase space are characterized by strong space time 
heterogeneities -the so-called dynamical hetcrogeneitics- 
which appear to be the by-product of an underlying first 
order dynamical transition in which the activity plays the 
role of an order parameter. A typical realization lies at 
a first-order transition point, which is thus characterized 
by the coexistence of competing families of trajectories. 
Active trajectories, in which a time and space extensive 
number of events where a particle escapes from its local 
cage, coexist with inactive trajectories in which local- 
ization of the particles within their local neighborhood 
dominates the dynamics. Without going further into the 
mathematical details of what the activity is, we present 
in figure [1] a snapshot of a the activity map in a three- 
dimensional KA mixture in three distinct situations. 

We now carry out our program and we begin by defin- 
ing an activity observable: given a set of particles inter- 
acting via a two-body potential V{vi —Tj), in which each 
particle i is subjected to a force Fj = — VV(ri — 
Tj), we introduce the observable: 
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The combination |F| + \ V ri • F, can be interpreted as 
the activity of a single particle. The quantity V e g in ([I]) 
appears in the study of Brownian interacting particles 
and measures the tendency for dynamical trajectories to 
evolve away from a given configuration. Indeed, it can 





FIG. 1: In these snapshots, the diameter of the particles quan- 
tifies the local activity: small points refer to mobile active 
particles, whereas large blobs refer to blocked inactive parti- 
cles. Red (blue) indicates that the activity is larger (smaller) 
than the median one. Left: Activity map in a typical config- 
uration from a trajectory with excess activity with respect to 
the average activity; there is a large number of mobile parti- 
cles though blocked particles also exist. Right: Activity map 
in a typical configuration from a trajectory with less activity 
than the average activity, large blobs dominate. 



be argued [T6| that the quantity exp(— pv c gdt) is pro- 
portional to the probability P(x, t + dt\x, t) that the sys- 
tem has stayed in its configuration between t and t + dt. 
This means that is the rate at which the system 

escapes its configuration. To understand how ([TJ is re- 
lated to this escape rate one writes P(x,t + dt\x,t) = 
(x\e~ \x) where H is the Fokker-Planck operator of 
evolution of the system. Using the detailed balance 
property of the dynamics, H can be symmetrized and 
this leads to P(x,t + dt\x,t) ~ sxp(— j3V e gdt), using 
standard path-integral representation of the propagator 
(x\e~ Hdt \x) (see [l6| for details). Besides, if one regu- 
larizes the Brownian dynamics in the form of hops on a 
lattice, our activity can be viewed as the continuum ana- 
log of the activity introduced in Lecomte et al. (l7j which 
counted the number of configuration changes undergone 
by a system over a given time interval. Perhaps more 
interestingly still, V c g also appears as the continuum 
limit analog of the dynamical complexity (a trajectory- 
dependent Kolmogorov-Sinai entropy §) which further 
clarifies its conceptual value. However, hand-waving ar- 
guments lead to a more concrete understanding of the 
value of the effective potential V e e- Indeed, in the ex- 
pression ([1]) minimizing the first term in the right hand 
side drives the system away from regions of phase space 
where forces are nonzero. In other words, it tends to favor 
mechanical equilibrium configurations (whether stable or 
metastablc). The second term in the right hand side of 
(fTj). can be positive or negative: if negative, it will favor 
local minima, all the more so as they are deep and steep; 
if positive, it will select local maxima. As a result, there 
will be two classes of trajectories, according to the val- 
ues of V e ff. It can be shown that the average equilibrium 
value of V e ff is (V e g) = —f3/4:J2i Pfi which is always neg- 
ative. For minimal (negative) values of V e ff, trajectories 
will explore deep energy basins, for maximal values of 
V e ff, trajectories will explore local maxima of the energy 
landscape. 
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FIG. 2: The dynamical free energy per particle ^^"^ is 
shown as a function of s for increasing sizes N = 45, 82, 
155, 250 and 393. While the curves collapse onto a single 
master curve at s < 0, a strong size-dependent behavior is 
observed at s > 0. On the s > side, the curves display 
a rapid variation at a typical value s c (N) which decreases 
as N increases, leading to the putative building up of a sin- 
gularity of %l>(s, N)/N at s = as N —¥ oo. The rescaled 
[ip(a, N) - ip(s c , N)]/N 1+a with a ~ 0.43 is shown in the in- 
set: on the inactive side s > the curves collapse onto a 
single curve, exemplifying the large N finite size scaling. Ap- 
plying the cloning algorithm requires the \s V e sdt\ product to 
be infinitesimally small (in principle), which at the left end of 
the graph and for the largest two values of N is hardly veri- 
fied (dashed lines). For these points, this product reaches .5 
which accounts for the increasing error. 

The total activity is defined as K(t) = f Q dt'V e s(t'). 
In terms of the local density p(x,£), our activity reads: 
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It is now apparent that T4ff involves three-body effective 
interactions, a fact that we a posteriori interpret by re- 
alizing that in order to escape from an energy basin via 
a collective rearrangement, multi-body interaction terms 
are needed. 

Given the activity we have just defined, we set out to 
determine the distribution of this quantity over a large 
number of realizations of the process, that is, given the 
dynamics is deterministic, by sampling over a large num- 
ber of initial states drawn from the equilibrium Boltz- 
mann distribution. Sampling the full distribution of the 
activity requires sufficient statistics beyond its typical 
values. For that reason we have preferred to work in 
an ensemble of trajectories in which the average activ- 
ity, rather than the activity, is fixed. In this canonical 
version, we consider a large number of systems evolving 
in parallel, and between t and t + dt we remove or add 
a system in a given configuration with a cloning factor 
equal to exp(— sV e s(t)dt). How to perform a numeri- 
cal simulation in a canonical ensemble of time realiza- 



tions was explained in Giardina et qZ.[15|'s cloning al- 
gorithm. Wc combined the cloning algorithm with the 
molecular dynamics simulation of the interacting parti- 
cles. Choosing a positive value of s allows one to focus 
on time realizations with a lower-than-average activity, 
while a negative s selects over-active trajectories, and 
s close to zero samples typical activity trajectories. To 
this day, no one has been able to endow the parame- 
ter conjugate to the activity s with a concrete physical 
meaning, so that the only physically realizable value of 
an s-ensemble is achieved for s — 0. As will shortly 
become clear, analyzing the s = properties gains enor- 
mously from being able to vary s away from 0, a theorist's 
privilege. The distribution of the activity K is fully en- 
coded in its generating function Z(s,t,N) = (e~ sK ), or, 
alternatively, in the corresponding cumulant generating 
function In Z. We define the intensive dynamical free 
energy tp(s,N) = lim^oo for any given finite size. 
From the knowledge of the function ip(s,N) we can re- 
construct the properties of the activity K. In numerical 
terms, in view of the cloning algorithm used, ip stands 
for the growth rate of the population of systems evolving 
in parallel. The details of the simulations arc as follows. 
We used the A-B mixture of [Io[ for samples of different 
sizes. The samples all had the approximate ratio 80/20 
for A/B particles. They were first prepared at equilib- 
rium at temperature T = 0.8 by coupling to a stochastic 
heat bath; the time step was dt = 2.10~ 2 . During a very 
long simulation at equilibrium, the N c clones where pre- 
pared. Then the cloning algorithm of Giardina et al. [l5| 
was performed with a small coupling to the heat bath, 
necessary to provide some stochasticity to the exploration 
of trajectory space, as explained in [3, [l8[ for determin- 
istic systems. In our simulations, the convergence was 
checked by varying the number of clones and the dura- 
tion of the simulations. For N = 45,82 the best results 
were found for N c = 1000 and r = 10r a (r Q is the relax- 
ation time); for N = 155,250,393 the best results were 
found for N c = 500 and r = 20r Q . Fluctuations in the 
space of clones would have required for the largest sys- 
tem sizes that we substantially increase N c as the tails of 
the distribution are sampled (extreme values of s). 

We now present the central result of the work, namely 
the dynamical free energy ip(s,N) as a function of s for 
increasing system sizes N. Figure [5] shows that the dy- 
namical free energy per particle ^^ jv - ) settles to a well- 
defined limiting value for s < as N increases. This 
means that for s < there exists a thermodynamic limit. 
On the s < side, where the activity is above its average 
value at s = 0, which we interpret in terms of active tra- 
jectories, the physical properties are not far from those 
at s — as predicted from the Gibbs equilibrium dis- 
tribution. However, a dramatic change in behavior is 
observed for positive values of s. Not only do the dy- 
namical free energies exhibit a strong size dependence, 
but a severe change in behavior is also observed on that 
side of the parameter space. For s > 0, the free energy 
rapidly increases, and the activity rapidly decreases. The 
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FIG. 3: The average activity per particle {V c a)/N is mea- 
sured in an ensemble of trajectories biased by a weight e~ aK 
and is plotted as a function of s. While the a < regime 
shows convergence to a size independent limit of an order 
comparable to its equilibrium value, the s > side confirms 
a strongly size-dependent behavior as V e fi abruptly decreases 
to negative values. This sudden drop in V e g indicates that 
the system freezes into inactive states. 
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FIG. 4: s c (JV) as a function of N. 

location s c (N) at which this rapid increase sets shows a 
strong A-dependence: s c (N) seems to decrease to zero 
as N increases. The precise location of s c (N) has been 
determined thanks to a spline method and by maximiz- 
ing tjj"(s); this is illustrated in figure 01 One notes that 
given the existence of error bars, there is an overall de- 
crease of s c (N) as N increases, but one cannot exclude 
that s c (N) converges to a very small positive value. As in 
other examples of first order dynamical transitions [BJ , we 
remark that s c {N) remains positive for all finite N. In- 
deed, ip'(0) = — (V e s) is proportional to N, which means 
that tp( s ) a l so scales as N in the vicinity of s = 0, and 
since s c (N) marks the transition to a regime where ijj(s) 
scales differently with N, one has s c (N) > 0. 

The interpretation in terms of active vs. inactive tra- 
jectories is easier when plotting the average activity as 
a function of s. This is done in figure [31 where the av- 
erage activity is measured independently: we see that in 
the vicinity of s c (N), the behavior of (V e s)/N changes 
abruptly from a smooth range of values at s < to 
a strongly negative range for s > 0. In this inactive 
regime, through rescaling (see the inset on figure [5]), we 



see that the free energy behaves in N 1+a with a ~ 0.43, 
in contrast to the purely extensive behavior of the active 
regime. A similar difference of scaling exponents in N is 
observed in KCMs in two dimensions @ at s > 0: ip(s) 
is of order L for the Fredrickson-Anderson model (FA) 
while of order 1 for the triangular lattice gas (TLG). This 
behavior is related to the geometry of the remaining ac- 
tive sites in the system: those divide up along the border 
of fully inactive domains in the FA model, while they are 
isolated in the TLG. Our finding that a is non- integer in- 
dicates that the KA mixture adopts configurations with 
non-trivial geometrical features for the inactive particles 
in the inactive phase due to the effective long-range in- 
teractions that develop in the s > s c states. Since large 
values of \V c s\ correspond to inactive histories, it is ex- 
pected that a > 0, as observed. We note that such finite- 
size behaviour of the activity with non-trivial exponents 
is known to occur in lattice KCMs [l9| , where the value of 
those exponents are directly related to the nature of con- 
figurations appearing in the inactive state (which differ 
from those of the active state). 

This allows us to claim that there exists a phase tran- 
sition from inactive to active states as s is varied. Our 
closely related goal was to identify the location of the 
transition. We have shown that s c (N) displays a strong 
size-dependence with an overall decrease for the range of 
N values explored. The two different well-defined sear- 
ings ~ _/V 1+Q (rcsp. <~ N) for s > (rcsp. s < 0) 
indicate that we have reached the large N asymptotics 
on our range of system sizes (see figure [2]) . We therefore 
conclude that the most likely value at which the transi- 
tion takes place in an infinite system is either s c (oo) = 
or a very small positive value close to 0. We were not 
able to perform simulations for larger values of N due to 
the very large computation times needed. 

The first order dynamic transition scenario observed in 
KCM's is thus confirmed in the atomistic model we have 
studied. As we mentioned earlier, the location of the 
phase transition along the s-variable axis is an extremely 
relevant issue since only the value s — is experimentally 
accessible. Note however that identifying a transition 
point at s = is possible only on the condition that s 
is varied across 0. Therefore, finding evidence for the 
transition occurring at s — renders an experimental 
observation of the transition a credible achievement. The 
lesson to be drawn from the present extensive simulation 
series is that, as expected, it is desirable to work on as 
small systems as possible, yet large enough to allow for 
collective effects in trajectory space to develop. 
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